Bogoliubov theory of interacting bosons on a lattice in a synthetic magnetic field 
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We consider theoretically the problem of an artificial gauge potential applied to a cold atomic 
system of interacting neutral bosons in a tight-binding optical lattice. Using the Bose-Hubbard 
model, we show that an effective magnetic field leads to superfluid phases with simultaneous spatial 
order, which we analyze using Bogoliubov theory. This gives a consistent expansion in terms of 
quantum and thermal fluctuations, in which the lowest order gives a Gross-Pitaevskii equation 
determining the condensate configuration. We apply an analysis based on the magnetic symmetry 
group to show how the spatial structure of this configuration depends on commensuration between 
the magnetic field and the lattice. Higher orders describe the quasiparticle excitations, whose 
spectrum combines the intricacy the Hofstadter butterfly with the characteristic features of the 
superfluid phase. We use the depletion of the condensate to determine the range of validity of our 
approximations and also to find an estimate for the onset of the Mott insulator phase. Our theory 
provides concrete experimental predictions, for both time-of-flight imagery and Bragg spectroscopy. 

PACS numbers: 



I. INTRODUCTION 

One of the main goals of experiments with cold atoms 
has been to extend to new contexts physical effects that 
are familiar from the study of condensed matter. From 
this point of view, a particularly important area has been 
the study of vortices in superfluid systems, previously 
considered in the context of type-II superconductors^ and 
superfluid helium^ When a trapped superfluid comprised 
of weakly interacting bosonic atoms is caused to rotate, 
it forms a lattice of vortices, a remarkably clear demon- 
stration of the quantization of circulation!^ 

In a corotating reference frame, the (neutral) atoms ex- 
perience a Coriolis force of the same form as the Lorentz 
force on charged particles in a magnetic fields This mag- 
netic analogy has inspired considerable experimental ef- 
fort to achieve effective magnetic fields large enough to 
reach the quantum Hall regime, and corresponding the- 
oretical work to extend the theory of the quantum Hall 
effect to the context of trapped bosons^ 

Much recent effort has been directed towards combin- 
ing effective magnetic fields and optical lattices, both to 
increase the stability of the experimental systems, and 
because of interesting physical effects expected in the 
presence of a lattice. The first experiments with rotating 
lattices used masks to produce parallel beams, whose sub- 
sequent interference formed the optical lattice potential. 
Mechanical rotation of the mask caused the interference 
pattern to rotate, and a density of vortices comparable to 
the density of lattice sites was achieved^ Mechanical in- 
stabilities limited the lattice strength, however, restrict- 
ing to the regime where the vortices are weakly pinned^ 

More recent experiments^ have replaced the mask with 
an acousto-optic modulator, allowing for considerably 
deeper lattices and lower temperatures. For a sufficiently 
deep lattice, a single-band approximation becomes valid, 
and the Bose-Hubbard mode&£ gives a good description 
of the physics. Experiments with Raman lasers^— allow 



for effective magnetic fields without rotation, instead us- 
ing coupling between internal atomic states to imprint 
the required geometric phases. A static optical lattice 
can be applied to such an arrangement, giving an effec- 
tive magnetic field within the laboratory frame. Various 
other proposals have been made to induce phases directly 
within the lattice*^— some of which have the advantage 
of producing a perfectly commensurate flux density. 

A major goal of these experiments is, as noted above, 
to reach the quantum Hall regime yi2r— in which the ef- 
fective flux density (in units of the flux quantum cf>o = 
2irh/Q where Q is the effective particle charge) is com- 
parable to particle density. Continuum systems in this 
regime exhibit a sequence of incompressible phases, the 
integer and fractional quantum Hall states, and, in the 
presence of a lattice, closely related physics has been pre- 
dicted in certain areas of the phase diagram^ 

Systems of bosons also support compressible superflu- 
ids, which are more closely related to the phases in the 
absence of a magnetic field, and are likely competitors 
with the quantum Hall states in the phase diagram. De- 
tailed study of these phases is beneficial for the search 
for quantum Hall physics, both to calibrate experiments 
and to understand this competition, but their nontrivial 
properties and phenomena, especially in comparison with 
conventional superfluids, make them worthy of consider- 
able interest in their own right. 

These properties are inherited from the remarkable 
structure of the corresponding noninteracting problem, 
a single particle moving on a tight-binding lattice in the 
presence of a uniform magnetic fielder— The density 
of states exhibits a fractal structure known as the 'Hofs- 
tadter butterfly' (see Section Hi CI and in particular Fig- 
ure^, with the spectrum depending sensitively on a, 
the magnetic flux per plaquette of the lattice (measured 
in units of 4>q). For rational a = p/q (with p and q co- 
prime), there are q bands, and each state is g-fold degen- 
erate. This degeneracy can be understood in terms of the 
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'magnetic symmetry group'^ which takes into account 
the modification of the symmetries inherent in making a 
particular choice for the gauge potential. 

The present work introduces a theoretical approach 
that starts with the noninteracting Hofstadter spectrum 
and treats the superfluid using Bogoliubov theory. This 
gives a mean-field description of the superfluid phase of 
the Bose- Hubbard model, in which the interplay between 
the 'magnetic' vorticity and the lattice potential (analo- 
gous to pinning effects for a weak lattice^^) leads to real- 
space density modulations ('supersolidity'). The theory 
furthermore provides a consistent expansion in terms of 
both thermal and quantum fluctuations, and describes 
the quasiparticle spectrum above the condensate, which 
combines features of the Hofstadter butterfly and the 
Goldstone mode characteristic of superfluid order. We 
present detailed calculations within the superfluid phase, 
based on the microscopic Hamiltonian, which include pre- 
dictions for both time-of-flight and Bragg spectroscopy 
measurements. A brief outline of our methods and re- 
sults has been presented elsewhere^ 

Bogoliubov theory has previously been applied to un- 
frustrated superfluids in optical lattices^ and also more 
recently to a proposed 'staggered-flux' modelM- which 
shares some features of the present analysis. Duric and 
Lee 3 ^ have applied a spin-wave analysis to the system 
considered here, using a real-space perspective that pro- 
vides results consistent with ours. A related analysis has 
been applied to fermions^i which similarly show spatial 
order in the paired superfluid phase. 

The real-space structure of the condensate also bears 
many similarities to other systems with phase coherence 
in a magnetic field. In particular, superconducting lat- 
tices with applied fields support states with current pat- 
terns similar to those discussed in Section Mil below ^ 
and experiments on these systems have clearly demon- 
strated the effect of the Hofstadter spectrum on the su- 
perconducting transition^ Frustrated Josephson junc- 
tion arrays ; 34 ' 35 where the charging energy can become 
comparable to the Cooper-pair tunneling amplitude, pro- 
vide a close analogue of the present system, and imple- 
mentations using cold atoms have been proposedi 36 ' 37 

Previous work on this system has also considered the 
Mott insulator that should exist for strong interactions 
and commensurate density) 7 ^ ' 26 ' 38 ' 39 This phase, which 
is very similar to its analogue in the absence of a mag- 
netic field, is favored when the Hubbard U interaction 
suppresses number fluctuations and eliminates the co- 
herence between neighboring sites of lattice. The tran- 
sition between this phase and the superfluid has been 
studied using the Gutzwiller ansat a 26 ' 38 " — and effective 
field-theory methods^ with the main effect being an en- 
hancement of the insulator phase due to the frustration 4 ^ 
of the hopping. 

In this work, we restrict consideration to uniform mag- 
netic fields. By modulating the effective field in real 
space, it is possible to produce states with nontrivial 
topological properties^ 5 - These include topological in- 



sulators and metals, which display similar physics to 
the single-particle states underlying the phenomena de- 
scribed here. Recent work has also addressed systems of 
interacting bosons in spatially varying magnetic fields4£ 

Although the framework that we present has broad 
applicability, our specific model includes several simpli- 
fications. First, we use a single-band Hubbard model, 
incorporating only the lowest band of the optical-lattice 
potential. (One effect of the magnetic field is to split the 
tight-binding spectrum into multiple 'Hofstadter bands', 
and these are included exactly in our approach.) While 
this is certainly not a valid approximation for the shallow 
lattices of earlier experiments^ the single-band limit is 
approached by more recent work with rotating lattices^ 
and should be readily achievable in combination with 
Raman-induced gauge potentials . 10 ' 11 

We also assume a spatially uniform system with ex- 
actly rational a. The main effect of the finite trap size is 
to wash out the small-scale fractal structure of the Hof- 
stadter butterfly; 13 ' 24 so we expect most of our results 
to be valid for any a. The experimental consequences 
of nonuniformity and incommensurate magnetic flux will 
be addressed in more detail in future work. 



A. Outline 

Our approach is as follows: In Section|Hl we use a sym- 
metry analysis to find the full spectrum of single-particle 
states, described by bosonic annihilation operators a^ 7 , 
defined below in Eq. (|I8|) . We then apply the Bogoliubov 
ansatz; 47 ' 48 

clm-i = A tl (2ir) 2 8 2 (k) + a k e~ ( , (1) 

where is a set of c-numbers giving the condensate 
order parameter, and the operators dke-y describe the 
residual bosons outside the condensate. An expansion in 
terms of these operators is then developed, which in phys- 
ical terms is an expansion in fluctuations around mean- 
field theory. In Section [TTTT we address the lowest-order 
term, which involves only the constants Ag^ and leads 
to a time-independent Gross-Pitaevskii equation for the 
condensate configuration. 

Higher-order terms, describing the spectrum of low- 
energy excitations and their interactions, are treated in 
Section IIVI In this work, we truncate the expansion 
at quadratic order, leading to a theory of noninteract- 
ing Bogoliubov quasiparticles. For this approximation to 
be reasonable, one requires a sufficiently low density of 
quasiparticles, and we determine the regime of validity 
by calculating the condensate depletion. 

In Section [Vj we describe the consequences for experi- 
ments, giving predictions for both time-of-flight imagery 
and Bragg spectroscopy. We conclude with discussion in 
Section I VII Some details of the Bogoliubov transforma- 
tion and the symmetry analysis of the interacting spec- 
tra are given in Appendices [X] and [B] In Appendix [C] we 
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provide more details and analytic results for the simplest 
case, a = |. 

B. Hamiltonian 

The single-band Bose-Hubbard model in the presence 
of a static magnetic field can be written-^ as 

(ij) 3 

where bj and rij = bjbj are the annihilation and num- 
ber operators respectively for site j. We treat a two- 
dimensional square lattice, and assume hopping with 
amplitude t only between nearest neighbors, denoted 
(ij). The second term gives an on-site interaction with 
strength U ; the approach described here can straightfor- 
wardly be extended to incorporate more complicated in- 
teractions, including between particles on different sites. 
(Note that we treat electrically neutral particles in the 
presence of a synthetic magnetic field, and we will as- 
sume that there are no long-range interactions.) 

The phase on the directed link from site i to site 
j, to be denoted i— >j, can be expressed^ , in terms of the 
magnetic vector potential *4. as 

- Q I ' dr ■ A(r) . (3) 

J Xi 

The integral is to be taken along the straight-line path 
between the positions Xi and Xj of the two sites; one has 
3>ij = —Qji- While the individual phases depend on the 
choice of gauge, the lattice curl 




is gauge- independent, where the sum is over links enclos- 
ing a given plaquette □ in a counterclockwise sense. 

The integral in Eq. (j4]) is simply the magnetic flux 
through the plaquette, given by \J3\a 2 , assuming a uni- 
form magnetic field B = V x A. perpendicular to a square 
lattice with spacing a. The dimensionless flux per plaque- 
tte is defined by a = \B\a 2 /<f)Q, and completely specifies 
the effect of the magnetic field. We will henceforth use 
units where h = a = 1, and take the effective charge on 
the bosons as Q = +1, so the flux quantum is simply 
<j>o = 2-7T, and Eq. @ becomes 

ijOD 

It is convenient theoretically, and also for recent ex- 
periments with effective gauge potentials ) 10 ! 11 to use the 
Landau gauge, in which the vector potential is parallel 
to one of the square lattice axes. We take the vector po- 
tential along y, which is conventional in the theoretical 
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FIG. 1: An illustration of the vector potential &ij in the Lan- 
dau gauge, and the symmetries defined in Section [II Al The 
number of arrowheads on each link of the square lattice gives 
the value of &ij on the directed link i—tj, in units of 2na. 
The configuration is one possible choice obeying Eq. ([5]): go- 
ing around any square plaquette in a counterclockwise sense, 
the number of forward arrows minus backward arrows is +1. 
In this choice of gauge, ^ vanishes on all links in the x di- 
rection. Any choice of gauge reduces the full symmetry of 
the square lattice, and the translation and rotation opera- 
tions shown are only symmetries when accompanied by phase 
factors (see Section Til A} . Reflections must be combined with 
time-reversal operations to preserve the direction of the ap- 
plied magnetic field. 

literature, but it should noted that in the recent experi- 
ments of Lin et al. , 10 ' 11 it is instead aligned with x. Our 
choice of gauge is illustrated in Figure [T] 

With this choice, one has e 1 * 3 '-^* = 1 and e'*w+* = 
c 2max 5 ; where (for example) j + x denotes the site adja- 
cent to j in the +x direction. The kinetic term 7i t then 
becomes 

«* = -*E[ fc 5+« 6 i+ 6 J-* 6 i+ w " 6 5+* b i 

+ u- x ib]_-b j ], (6) 

where oj = e 27na . Note that the phases act to 'frustrate' 
the hopping, so that for noninteger a it is not possible to 
minimize the kinetic energy on every link of the lattice 
simultaneously^ 

The transformation to an alternative gauge is imple- 
mented by applying a spatially-varying phase rotation 
to bj. For example, in the symmetric gauge, which is 
more natural for a rotating lattice, one defines opera- 
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tors bj — w x i y i/ 2 bj, so that the gauge field becomes 
e i*jj+x — lj~ v ^ 2 and e'*"+ s = uj Xj / 2 . The number 
operator rij = bjbj is invariant under any such gauge 
transformation, as required for a physically observable 
quantity, and so the interaction energy Hu is also invari- 
ant. 

As noted above, we assume precisely rational a = p/q 
(with p and q coprime), so that u q — 1, and the unit cell 
of T-L is q x 1 sites. (In the symmetric gauge, the unit cell 
is considerably larger, containing 2q x 2q sites.) 

II. SINGLE-PARTICLE SPECTRUM 

We begin by describing in detail the spectrum of the 
single-particle kinetic term Ht, given in Eq. (J5]). While 
these result s 21 ' 23 ' 24 are well established (and have been 
for many decades), we will present them here in some 
detail, in order to introduce the concepts and formalism 
that will be central to our subsequent analysis. 

The general structure of the spectrum for any a — p/q 
can be determined by considering the "magnetic symme- 
try group" (MSG; also known as the "projective symme- 
try group"). 23 ' 49 This arises because any choice of gauge 
necessarily reduces the physical symmetry of the lattice 
and leads to a Hamiltonian that does not commute with 
the standard spatial symmetry operators. For exam- 
ple, while the (uniform) magnetic field is invariant un- 
der translation by a single lattice site in any direction, 
T~it is manifestly asymmetric under translations in the x 
direction. 

One can nonetheless define a group of operators in 
one-to-one correspondence with the physical symmetries, 
that obey the multiplication table apart from phase fac- 
tors, and that commute with Tit (and indeed the full 
Hamiltonian %). We define this group by specifying op- 
erators corresponding to the elementary translation, ro- 
tation, and reflection transformations from which the full 
group can be constructed. 

A. Magnetic symmetry group 

We first define the elementary translation operators l~ x 
and Tj,, illustrated in Figure [TJ in terms of their commu- 
tation with the annihilation operator bj . With our choice 
of the Landau gauge, the Hamiltonian is symmetric un- 
der translations in the y direction, and so one can define 
T y by 

'/•A. = h 3+y% ■ ( 7 ) 

In contrast, the x-dependent phase factors in Eq. (JB]) im- 
ply that a pure translation does not commute with "H t , 
and we instead define T x by 

T x bj = b j+x T x u- y * . (8) 

The gauge field has periodicity q in the x direction, so 
the combination T£ obeys T x q bj = bj+ qx T x q - 



The multiplication relations of the MSG are equal up 
to phase factors to those of the ordinary spatial group. 
With total particle number N = J^j n j i one finds 

T x Ty = T v T x oj N , (9) 

by induction, starting from the (totally symmetric) vac- 
uum state. While the phase factors associated with in- 
dividual operators are dependent on the choice of gauge, 
this relation is gauge independent. 

Besides translations, it is useful to consider rotation 
and reflection operations. We define the unitary operator 
1Z giving a rotation by 90° counterclockwise about the 
site at the origin (see Figure [J): 

Rb, 7C. (10) 

where j — > WLj under the rotation (xm_j = ~Vj, y«.j = 
The phase factor again ensures that [R., Tit] = 0. 

For reflection operators, the situation is somewhat dif- 
ferent, since the magnetic field explicitly breaks chirality 
and time-reversal symmetry. While reflection reverses 
chirality and so is not a symmetry of the Hamiltonian, a 
combination of reflection and time reversal restores the 
appropriate sense of circulation and remains a symme- 
try in the presence of the field. One can therefore de- 
fine antiunitary operators corresponding to such com- 
binations; we define T x and T y for the transformations 
obeying xi^j = —Xj and yi j = —yj respectively. No 
phase factors are required in these cases. 

It is also useful to define the combinations V = TZ 2 = 
I x X y and X xy = lZI y . The former gives inversion about 
the origin, x — > — x, which in two dimensions is a proper 
rotation and hence represented by a unitary operator. 
The antiunitary operator I xy gives reflection in the line 
y = x. 

Note that the interaction term Jiu is a function only of 
the gauge-invariant combination rij = bjbj and so is un- 
affected by the phase factors included in expressions such 
as Eq. ([S]). Any interaction term with the full symmetry 
of the lattice, such as the explicit example in Eq. @ , is 
therefore invariant under the magnetic symmetry group. 



B. Momentum-space operators 

While the unit cell of the Hamiltonian T-L t contains q 
sites and hence gives a reduced Brillouin zone, it is useful 
to construct momentum space operators based on the 
full lattice Brillouin zone Q3l: — 7t < k x ,k y < tt. The 
momentum-space annihilation operator bk is defined by 

'</, >> : " : ''V (ii) 

3 

so that the commutator is given by 

[b k ,bl] = (2n) 2 S 2 ([k-k%J. (12) 
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Here and throughout, we use the notation [k]<s to denote 
the momentum k reduced to the Brillouin zone 03 by the 
addition of an appropriate lattice vector. We also use the 



shorthand notation [x] q = x mod q. 

Written in momentum space, the Hamiltonian T-Lt is 



Hi 



fees 



d 2 fc 

WT 2 



2co S k x b\b k + c [k "b\ k+x] b k + c lk yb\ k _ x] b k 



(13) 



where X = lixax. The enlarged unit cell (ijx 1 sites) of the Hamiltonian allows mixing between momentum states 
that coincide when reduced to the magnetic Brillouin zone 03m: — tt < k y < tt, —tt/q < k x < Tt/q- 
The operators bk commute with T yi 



Tybk 



<b k %, 



but the phase factor in the definition of T x causes it to mix momenta, 

T x bk = e lfex 6[ fc+ y] 3 , L 7i , 



(14) 



(15) 



where Y = 2nay. Since T x commutes with H t , this implies degeneracies in the single-particle spectrum between 
points separated by Y. To make this transparent, it is convenient to define the doubly reduced Brillouin zone 03n: 
-n/q < k x , k y < ir/q. 

Any point within 03l can be specified as [k + IX + nY"]<8 L , where k £ 03n and n and I are integers from to q — 1, 
so %t can be rewritten as 



f d 2 k q ~ 1 q ~ 1 



(16) 



1=0 n,n'=0 



r 



where the q x q matrix H(fe) has elements (for q > 2) 



H n n' \k) 




if n' = n 
if n' = In 



if n' = [n — l] q 
otherwise. 

(17) 

(If q = 2, [n+l] g = [n— l] g and Hqi = Hiq = —2t cos k y .) 

Noting that H nn ,([k + £Y]<sJ = ^ n - n '^H nn ,{k), we 
diagonalize "H f by writing 
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ip~fn(k)a k i 7 , (18) 



where ake-y is the annihilation operator for a single- 
particle state labeled by band index 7 £ {1, . . . , q}. For 
each k £ 03 n, i ! ^n{k) is an eigenvector of H(fe) with 
eigenvalue e 7 (fc). The q bands can be understood from 
the 'folding' of the Brillouin zone due to the reduced 
translation symmetry of %. 

It should be noted that the annihilation operator a k ij 
has momentum k when referred to 03n, but momentum 
k + £Y in the magnetic Brillouin zone 03m- For each 
k £ Q3n, the eigenvectors ipjn(k) corresponding to dif- 
ferent bands are orthogonal, so the operators akt-y obey 
canonical commutation relations, 



1 k£y a k'£'~ 



,} = (27r) 2 S 2 ([k-k% N )5 u -S^. (19) 



The single-particle Hamiltonian can finally be rewrit- 
ten as 



Ht = 



fees. 



d 2 fc 

(2nf 



q-l 



kt~t a kl~t ■ 



(20) 



e=o 7 



The single-particle energy e 7 (fe) is independent of £, so 
every state is (at least) g-fold degenerate. The band la- 
bels 7 can be arranged so that e 7 (fe) < e 7 +i(fc) for every 
fc and e 7 (fe) is a continuous function of fc. 

Using Eqs. (fT4|) and (fl~5| . one finds the effect of the 
operators T y and T x as 



T* ik £ T' 

T y ak£y = e v u> ak£ 7 ly 
T x ak£~ ( = e lkx ak[e+i] q -yT x 



(21) 
(22) 



The operator T x therefore transforms one degenerate 
single-particle state into another, and it is this symmetry 
that enforces the degeneracy. 

Determining the effects of the rotation and reflection 
operators 1Z, I Xl and I y is somewhat more involved^ 
Considering first 72., its commutation with Ji t implies 
that IZcLkt-y'R^ can be written in terms of fl(Hfe)f 7 in the 
same band 7. With an appropriate choice of the arbi- 
trary phase of the eigenvectors ip in at the points k and 
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Rfc, one finds 

Tlake-, = — c -1 *"' ^2 w " a (Rfc)£'7^ 1 (23) 

v 

where the phase 7 is independent of fe. The requirement 
that 1Z 4 = 1 implies that </> 7 is a multiple of ir/2 for all 
bands 7. One similarly finds 

%cQ-kti = a(i yk )[^i] q7 l x (24) 
IyOfc£ 7 = e 2l * T a(i x k)e-fl y (25) 
Pa M7 = e- 2i ^a ( _ fe)M87 P. (26) 

Note that e _2l< ^ T = ±1 is real. Similar expressions, such 
as 

W(Rfc) = ■^e- i *'J]« flB V'( fc ) - ( 27 ) 

relate the eigenvectors ipy n (k) at symmetry-equivalent 
momenta fe. 

Applying these operators to the Hamiltonian in the 
form of Eq. (f2"0")) immediately shows that the single- 
particle dispersion e 7 (fc) is symmetric under the corre- 
sponding transformations of the momentum fe. For ex- 
ample, e 7 (Rfc) = e 7 (fe), which implies that the dispersion 
has the full four-fold rotation symmetry of the lattice, de- 
spite the reduced symmetry of %t ■ 



C. Spectrum 

In summary, for a = p/q, the single-particle spectrum 
consists of q bands, labeled by 7, with each state g-fold 
degenerate. These degenerate states have energy e 7 (fc) 
and momentum [fe + £Y]<g M . with I £ {0, ...,q — 1}. 
Figure [5] shows the 'Hofstadter butterfly'^ a plot of the 
allowed single-particle energies e 7 (for any momentum 
fe) as a function of the flux a. The plot has a fractal 
structure 2 ^ that is sensitively dependent on a, and for 
clarity only rational a = p/q with q < 10 have been 
included. For each a = p/q, points mark the top and 
bottom of each of the q bands. 

As can be seen in Figure [2 most of the bands are 
separated by nonzero gaps. The only exceptions are the 
two central bands for q even, which touch exactly at the 
point of zero energy. These occur at fe = for q an 
integer multiple of 4, and at the corner of *8n, fex = 
— x + ^y, otherwise. For both, the spectrum has a linear 
'Dirac-cone' dispersion near the degeneracy points 

In all other cases, including the lowest band for all a, 
the dispersion is quadratic near its minimum. The mini- 
mum of the lowest band always occurs at fe = 0, as can 
be shown using the Perron-Frobenius theorem, and the 
effective mass near this point can be found by perturba- 
tion theory for small fe. In all cases, the coefficients are 
equal in the x and y directions (i.e., the effective mass 
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FIG. 2: The Hofstadter butterfly,— a plot of the single- 
particle energies e 7 as a function of flux 01. The butterfly 
has a fractal structure, but only a finite set of a = p/q can 
be plotted; for clarity, we restrict to q < 10. Points mark the 
top and bottom of each of the q bands, which become increas- 
ingly narrow as q becomes larger. In the limit q — > 00 with p 
fixed, or equivalently a < 1, the low-lying bands become the 
Landau levels of the continuum. 

tensor is proportional to the unit matrix), a straightfor- 
ward consequence of the symmetry 1Z. Figure[3]shows the 
lowest band of the noninteracting dispersion, ei(fe), for 
a = I and jk In both cases, the dispersion is quadratic 
and isotropic near the top and bottom of the band. 

D. Interactions 

The operators defined in Section Hi Bl are chosen to 
diagonalize the kinetic energy operator "H t . To incorpo- 
rate the effects of the interaction term Hjj, this must also 
be expressed in terms of these operators. This involves 
the straightforward process of substituting Eqs. (fTTj) and 
(|18p into Hu, and can be performed for any choice of 
interaction. Our explicit calculations are for the on-site 
Hubbard interaction in Eq. ^ . appropriate to bosons in 
a deep optical lattice. 

A general quartic interaction can be written in terms 
of the operators aui^ as 

= ^2 X! Ua L^l7i a L«272 a fc 3 «373 a fc4«474 ' 

Jfc l'" fc -1 ly-U 71-74 

(28) 

where the coefficient u is a function of the four sets of 
indices fe, I, and 7. It can be chosen symmetric under 
exchange of the first two (1 2) or last two sets (3 H 4), 
and the requirement that T-Ljj be hermitian implies that 
u(3,4 ) l,2) = u*(l,2,3,4). 

A translation-invariant interaction conserves momen- 
tum, so the coupling coefficient u is nonzero only if 

[(k 1 + k 2 -k 3 -k4) + (h+£ 2 -£ 3 -h)Y} <3 =0. (29) 
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Note that this allows for umklapp processes where the net 
momentum is zero only when reduced to Q3n- Factoring 
out (27r) 2 (5 2 ([fci + fc 2 - fc 3 - fc4]<B N ) S ives S: 



2.91 



-2.92 
-2.93 



-2.95 



2.9(i 



FIG. 3: Contour plot of the lowest band of the noninteracting 
single-particle spectrum je 1 (k), for a = - = ^ (left) and ^ 
(right). In both cases, the spectrum is plotted in the magnetic 
Brillouin zone 25m, in which — — < k x < — and — n < k y < n. 
Points separated by momentum Y — 2nay are degenerate; 
identifying these gives the doubly reduced Brillouin zone 23n, 
— ^ < k x ,k y < indicated by the horizontal dashed lines. 
Note that the lowest band is considerably narrower in the case 
a = i, as is also evident in Figure [2] 




u({k},{£},{j}) = -6 W E ^{n}^ 1 n 1 (fel)^n 2 (fe2)^ 73n3 (fe3)V' 74 „ 4 (fe4)w" 1&+ "^-^ 3 -^, (30) 

ni---ft4 

where Sin and 5i n y denote Kronecker deltas enforcing Eq. (|29|) and a similar constraint on The complicated 

structure of the single-particle states gives u a nontrivial dependence on the momenta k of the interacting particles, 
despite the choice of a purely on-site interaction. 

Note that u is, up to a phase factor that is only nontrivial for umklapp processes, only dependent on two of the £'s: 

u{h, £ 2 , 4, 4) = e - i ^( fel + fe2 - fe 3-^)-*u(0, \h - £x] ? , \h - £i] q , [h - k] q ) , (31) 

where p is the modulo-q reciprocal of p, the integer such that [pp] q = 1 and < p < q. This identity, and the 
momentum-conservation constraint of Eq. (|29[) are consequences of the symmetry of the interactions under T x and 
Ty Further constraints on the coefficients u result from the symmetry properties of the eigenvectors ip in (k) under 
rotations and reflections. For example, requiring that V commutes with Ti.\j and using Eq. (|26|) gives 

e -2i £* =1 ^ u ( fel . . . fe4) 4 . . . i 4) 7l . . . 74 ) = . . . _ fe 4j [_4] g . . . 7l . . . 74 ) . (32) 

I 



This implies that u is odd in momentum, and hence van- 
ishes for fei.„4 = 0, for certain combinations of £i...4 and 
71...4. A detailed discussion of the restrictions imposed 
by symmetries has been presented by Balents et al^ in 
a different context. 



III. MEAN-FIELD THEORY 



Having described the spectrum of noninteracting par- 
ticles, we now turn to the effects of interactions on the 
many-body physics. As described in Section ITAl our ap- 
proach will be based on the Bogoliubov theory, using the 
ansatz Eq. ([T]) and performing an expansion in powers 
of the fluctuation operators. The Bogoliubov ansatz can 
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be viewed as a statement about correlation functions in 
the superfluid phase, with the first term in Eq. ([1]) giv- 
ing the one-point correlation function, the condensate or- 
der parameter (ofe^ 7 ) = A^ 7 (27r) 2 <5 2 (A;). The zeroth-order 
term in the Bogoliubov expansion is given by neglecting 
higher-order connected correlation functions. 

This mean-field theory can be viewed as a special 
case of that derived by using the Gutzwiller ansatz, 
which assumes a state FT ■ \ipj) that is factorizable in real 
spacei 26 i 38 ~ — In general, one allows 1^) to be an arbi- 
trary state within the on-site manifold, but our ansatz 
assumes a bosonic coherent state and is appropriate only 
within the superfluid phase. 

Substituting the mean-field ansatz into the Hamilto- 
nian, written as in Eqs. (T2TJ]) and (|28p. gives the energy 
density 

+ u({0},{t},{l})At ni Al l2 A e3l3 A eili , (33) 

m,{7} 

where the interaction strength u is evaluated with all 
momenta equal to zero. [This expression has been di- 
vided by a factor of (2tt) 2 S 2 (0), corresponding physically 
to the system volume.] The mean-field condensate con- 
figuration can be found by minimizing ho with respect 
to Agy. The resulting equation can be viewed as a time- 
independent Gross-Pitaevskii equation for the conden- 
sate wavefunction in momentum space. 

The corresponding real-space wavefunction can be 
found using Eqs. (TTTj) and (JT5J) and is given by 

(b,) = J2 u™^~ nt J2 VV(0)^ 7 • (34) 

In 7 

This is in general a function of [x] q and [y] q , and so gives 
a q x q site unit cell in real space. To this order, the par- 
ticle density is simply given by (rij) = \{bj)\ 2 . Note that 
the presence of A^ for nonzero I implies that the con- 
densate contains components for fc = IY ^ and hence 
that spatial symmetry is broken. Within this mean-field 
theory, this spatial order develops simultaneously with 
the breaking of phase-rotation symmetry, and is a sim- 
ple consequence of the degeneracy in i. (In fact, the 
finite-temperature transition in two dimensions is of the 
Berezinskii-Kosterlitz-Thouless type, and so this partic- 
ular result is not necessarily reliable.) 

The configuration of currents within the superfluid 
phase can be calculated using the gauge-invariant cur- 
rent operator for the link 

Jij = itj* i *b\b i +h.c. (35) 

Figure [4] shows the currents ( Jij ) in the mean- field con- 
densate configurations for a = i with 2 < q < 5; they 
have the same q x q unit cell as the condensate wave- 
function. For the larger values of q, particularly a — -|, 



-• 




FIG. 4: Example mean-field condensate configurations for 
a — i (top left), | (top right), i (bottom left), and i (bot- 
tom right). The blue arrows show the direction of the current 
{Jij) on each link i— >j of the lattice, and their lengths indi- 
cate the magnitude (the length scale is not consistent between 
different values of a). The black points show the positions of 
the lattice sites i and have area proportional to the density 
(rii). In each case, the condensate reduces the spatial symme- 
try of the square lattice, and is one member of a discrete set 
of degenerate configurations related by the action of the bro- 
ken symmetries. The degeneracies and residual symmetries 
are listed in Table U The quantitative details, but not the 
symmetries, depend on the interaction strength U ; the plots 
show the case U <(. 



these resemble Abrikosov lattices -X the plaquettes with 
low density and high current circulation can be viewed 
as containing vortices. The symmetry properties of these 
configurations are listed in Table |U 

An analysis of the patterns that are allowed for general 
interactions and various values of q has been given by 
Balents et aLj^S who considered the same problem in a 
different context. In the present case, it is valid to assume 
purely on-site interactions, allowing the ordered states to 
be determined unambiguously. 

Minimization of ho with respect to Ag 1 is equivalent 
to minimizing with respect to the real-space condensate 
wavefunction or vortex configuration. The latter perspec- 
tive is more appropriate in the continuum, whereas here 
the lattice potential provides a strong pinning potential 
that simplifies the momentum-space approach. 

Our ansatz for the condensate configuration, Eq. (JTJ, 
also involves bands 7 other than the lowest. Occupa- 
tion of higher bands costs kinetic energy, increasing the 
first term of Eq. f|33|) . and so is disfavored when inter- 
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degeneracy symmetries 
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6 

16 
10 



T T V 2 , 

fyfx: 'V , -L-xy 



/v i 1 y 1 x -L-y 



iy I x i 
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TABLE I: Properties of the mean-field condensate configu- 
rations shown in Figure [4] including their degeneracies and 
unbroken symmetries. The vectors in the column labeled An 
are configurations minimizing ho in the limit of weak inter- 
actions, U /t — > 0, where the condensate is restricted to the 
lowest band, 7 = 1. The residual symmetry group is given by 
products of powers of the operators listed, along with T x and 
T y q , which are always preserved by the ansatz of Eq. JTJ). The 
symmetries T x , 7~ y , and R are illustrated in Figure [T] and the 
combination X xy = T X 1Z gives a reflection in the line y — x. 



actions are very weak. For stronger interactions, the en- 
ergy is reduced by smoothing out density fluctuations, 
which requires incorporating higher bands into the con- 
densate. This competition between kinetic and potential 
energy also allows for first-order transitions between dif- 
ferent local minima of ho as the interaction strength or 
mean density varies. We have not found any examples 
for q < 5, however, and their observation in experiments 
would anyway likely require considerable enhancements 
in stability and cooling. 

The mean-field energy hg is symmetric under the same 
transformations of Ag 7 as the full Hamiltonian is un- 
der transformations of aot-y, as discussed in more de- 
tail in Appendix [Bj As noted above, certain symme- 
tries are spontaneously broken by the condensate con- 
figuration, and the corresponding operators transform a 
given into a symmetry-equivalent degenerate con- 
figuration. The degeneracies of the patterns shown in 
Figure [4] are listed in Table HI The number of degener- 
ate configurations is in every case a multiple of q, as we 
prove in Appendix [B] The degeneracy in the ordering 
patterns allows for the possibility of real-space domain 
formation, which would not affect time-of-fiight images 
and would likely require more sophisticated in situ probes 
to confirmj 52 i 53 

It should be noted that the ansatz of Eq. ([1} implicitly 
excludes ordered states with larger unit cells than q x q 
sites. (Previous work using a real-space approach^ has 
suggested that this may happen for a = j.) It is straight- 
forward to include such states (with larger but finite unit 
cells) at the mean-field level, by allowing nonzero conden- 
sate amplitude at a discrete set of momenta [fe]gs N 7^ 0. 
This complicates somewhat the analysis that follows, and 



we will not treat this possibility further. 

The condensate configuration A^ also determines the 
occupation numbers in momentum space, and so can 
be used to predict the result of a time-of-flight expan- 
sion measurement, as discussed in detail in Section [V Al 
Briefly, the terms in the Hamiltonian Eq. (|13p mixing 
momenta differing by X imply Bragg peaks at points 
corresponding to momenta nX, while the nonzero con- 
densate amplitude A^ 7 for i ^ gives further peaks at 
IY . (It should be recalled that our axes are reversed 
from those of Lin et ah 10 ! 11 ) 



IV. BOGOLIUBOV THEORY 

The mean-field theory of Section Mil results from us- 
ing the Bogoliubov ansatz of Eq. (JTJ) and keeping only the 
lowest-order term in an expansion in terms of the fluctua- 
tion operators a&f 7 . To improve upon this theory and de- 
termine the spectrum for single-particle excitations above 
the condensate, we consider in this section the following 
order in the expansion. The terms containing a single 
operator vanish when the mean-field energy density ho is 
minimized, and so we next treat the quadratic terms. 

The quadratic part of the Hamiltonian can be 

conveniently expressed in matrix form, by combining cre- 
ation and annihilation operators into a column vector, 



(36) 



-(-k)e 7/ 



One can then write T-L^> as 



(2) 

where % c is a term that contains no operators and 
comes from a commutator. 

This expression can straightforwardly be generalized 
to allow for other choices of single-particle basis, by re- 
placing £ and 7 by a single generic index A. The matrix 
M(fc) can then be written as 

Mav (fc) = la [ex x> (*0 - ^5xx> ] + B xx , (fc) , (38) 

where exx' , the generalization of e Sa>, is not diago- 
nal in the general case, and 



d 2 fc 



£',77' 



+ nf\ (37) 



j 



B A v(fc) 



/ Au(k,0,k,0;X,\ 1 ,X',X 2 )A* Xi A X2 2u(k, fc, 0, 0; A, A', Ai, X 2 )Ax 1 Ax 2 
^ \ 2u*(k, -fc, 0, 0; A', A, X u X 2 )A* Xi A\ 2 4u(-k, 0, -fc, 0; A', Ai, A, X 2 )A* Xi A X2 



Ai A 



(39) 
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Similarly to ho in Eq. (pJ5f . M(fe) has contributions 
from both the kinetic and potential energy. The latter 
can be viewed as self-energy terms for the quasiparticles 
due to scattering with bosons in the condensate. They in- 
clude 'anomalous' processes in which a pair of condensed 
particles scatter from each other into an excited state and 
the reverse process where they return to the condensate. 
These result in the off-diagonal elements in Eq. (|39)) . 
giving terms in 'US 2 ^ that do not conserve the number 
of dkij quanta i 47 i 48 The standard Bogoliubov theory for 
zero magnetic field is recovered by taking q = 1, in which 
case the t and 7 indices are redundant. 

It is useful to consider M^^/y (fe) as a 2q 2 x 2q 2 ma- 
trix (for each k). Using the properties of u given after 
Eq. (|28p . one can show that this matrix is hermitian. 
The zero-momentum limit M(0) is the Hessian of the 
mean-field term ho (with respect to variations in and 
its conjugate) and so is a nonnegative-definite matrix. 
The single vanishing eigenvalue corresponds to the bro- 
ken U(l) symmetry of ho- For nonzero k, all eigenvalues 
of M(fe) are strictly positive. 



A. Bogoliubov quasiparticles 

To find the spectrum of quasiparticles, one must define 
a new set of annihilation and creation operators in terms 
of which 'HP^ is diagonal. Momenta (referred to 2$n) are 
not mixed in Eq. (f37|) . so the new operators are labeled 
by k, but since the condensate breaks symmetry under 
T y , £ is no longer a good quantum number. We there- 
fore define annihilation operators for these modes as , 
where £ £ {1,2,..., q 2 }. (In certain cases, there are un- 
broken translation symmetries, as shown in Table [I] The 
states can then be labeled by the eigenvalues of the cor- 
responding operators, as discussed in Appendix [Bj) 

In order to preserve the bosonic commutation rela- 
tions, [dfe£,c4./£/] = (27r) 2 <5 2 (fc — k')5^/, the transforma- 
tion between and dk^ must be symplectic^ 4 - details 
are given in Appendix [Al In terms of the new operators, 
the quadratic part of the Hamiltonian is given by 



d 2 k 



(2tt 



6=c 



(2ir) 2 5 2 



where 



I k£y 



- r fcf 7 «(_fc)f 7 



(40) 



(41) 



To this order, the system is therefore described by non- 
interacting Bogoliubov quasiparticles with annihilation 
operators d^ and energies £fc£ > 0. Because of the off- 
diagonal elements in M(fe), the quasiparticles are super- 
positions of particles and holes, with v< " 
spectively giving these components. 



^kl-< aild *M 7 r 



r, 
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FIG. 5: Quasiparticle dispersion (solid lines) and noninter- 
acting single-particle dispersion (dashed), both in units of 
hopping t, for a = - = i. An analytic expression for the 
spectrum for this case is given in Eq. (|C7|) of Appendix [C] 
The dispersions are plotted along a path in the reduced Bril- 
louin zone 5?n shown in the left inset. In the interacting case, 
U = 4t, the mean density is p = 1, and the real-space con- 
figuration is as shown in the right inset (see also Figure [4|. 
In both cases there are q 2 = 4 modes, including, in the in- 
teracting case, one Goldstone mode with linear dispersion. 
For C/ = 0, the modes are g-fold degenerate and have been 
shifted vertically by an arbitrary choice of chemical potential. 
At k = + the corner X of 23n, the modes meet at 
a point, with a linear dispersion. Such a 'Dirac cone' occurs 
whenever q is even. Other notable features of the interacting 
spectrum include a twofold degeneracy along the line from M 
to X, and the unshifted modes (relative to the noninteracting 
dispersion) from X to F. The former can be understood as a 
Kramers degeneracy due to the antiunitary symmetry under 
Txly, as discussed in Appendix [B] 



For vanishing interactions, the second term in Eq. (|38p 
is absent and the quasiparticle spectrum is identical 
to the single-particle spectrum e 7 (fc) described in Sec- 
tion III CI With nonzero interactions, the q-fold degen- 
eracy within each band is split and, for generic k, the 
spectrum consists of q 2 distinct modes. The quasiparticle 
dispersion for a = \ and ^ are shown in Figures [S] and [S] 
respectively, along with the noninteracting single-particle 
spectrum. (For clarity, only 6 out of the q 2 = 9 modes 
are shown in the latter case.) For general q, the diago- 
nalization of M(fc) must be performed numerically, but 
the simplest case is analytically tractable, and is treated 
in detail in Appendix [Cl 

As the figures show, the lowest energy approaches zero 
in the limit k — > 0, giving the Goldstone mode that 
results from broken U(l) symmetry. For small |fe| this 
'phonon' has a linear dispersion, and can be described in 
terms of long-wavelength fluctuations of the condensate 
phase. A low-energy theory of the mode can be found 
by allowing gradual deviations from the mean-field value 
of the phase and its conjugate density, as described in 
Section IC2l 
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FIG. 6: Quasiparticle dispersion (solid lines) and noninter- 
acting single-particle dispersion (dashed), for a = 2 = -. In 
the interacting case, U = 2t, the mean density is p = 1. In 
both cases, there are q 2 = 9 modes, of which only the lowest 
6 are shown. The dispersions are plotted along a path in the 
reduced Brillouin zone QJn shown in the left inset. Because 
the condensate configuration breaks the symmetry under 1Z, 
rotation by £, the quasiparticle dispersions along the lines 
from T to X and from F to X' are different. The interact- 
ing spectrum includes several (unavoided) level crossings, for 
example in the lower band near X', a result of the unbroken 
translation symmetry T y T x , as discussed in Appendix [B] 



The phase velocity c of the Goldstone mode, which 
is expressed in terms of the spectrum at k = in Sec- 
tion IA 3( is in many cases independent of direction, in- 
cluding both a = 5 and ^. This isotropy is a straight- 
forward result of the symmetry in both cases under 
X xy = IZly, as noted in Table HI This property, and 
other features of the spectra shown in Figures [S] and [51 
are discussed in Appendix [Bl 

The second term in Eq. (|4T)|) includes ?i c ' from 
Eq. (|37|) and represents the change in the zero-point en- 
ergy associated with the superfluid state. This term, 
coming from quantum fluctuations, is accompanied at 
nonzero temperature T by a contribution from thermally 
excited Bogoliubov quasiparticles, leading to a free en- 
ergy per site of 



A/ = 



fee®] 



d 2 fc x - 

(2tt) 2 ^ 



Tlog(l-e-« fc «/ T ) 



(42) 

where we use units such that fee = 1 . (Within the mean- 
field theory of Section [TOl all particles are in the conden- 
sate, so the entropy vanishes and the free energy is given 
by ho.) These contributions in principle allow a configu- 
ration with a higher mean-field energy ho to be selected 
because of its enhanced fluctuations and hence lower free 
energy ho + A/. It should be noted, however, that the 
degeneracy of the symmetry-equivalent condensate con- 
figurations discussed in Section IHII cannot be lifted by 
A/. 



The calculated spectra lead to important experimental 
predictions, as discussed below in Section [V] Occupa- 
tion of the quasiparticle modes, due to both thermal and 
quantum fluctuations, gives the structure of time-of-flight 
images away from the Bragg peaks mentioned previously, 
and spectroscopic methods should be able to measure the 
mode dispersions £fc^ directly (see Section IV Bp . 



B. Condensate depletion 

The ansatz of Eq. (JTJ and the expansion in powers 
of operators is in principle exact, with the higher-order 
terms leading to interactions between the Bogoliubov 
quasiparticles. Here, the series is truncated at quadratic 
order, an approximation that is valid provided that the 
quasiparticles remain at sufficiently low density for their 
interactions to be neglected. 

This criterion can be quantified by calculating the de- 
pletion of the condensate, equal to the quasiparticle con- 
tribution to the total particle number. This is found by 
expressing the number operator rij is terms of the quasi- 
particle operators dk(, summing over sites j, and taking 
the ensemble average. The mean particle density is then 
given by 



r( |2 



+ 1^111 + ^(^)1. (43) 



where tib(£) = (e^ T — 1) _1 is the Bose-Einstein distribu- 
tion function. The first term in Eq. (|4"3"|) is the condensate 
density, and is simply the spatial average of the mean- 
field density calculated in Section IIII1 while the second 
term gives the average density of particles outside the 
condensate. The relative magnitude of these two terms 
gives a measure of the significance of fluctuations, and we 
use the ratio of the second to the first as our definition 
of the depletion. 

Figure shows the depletion for function 
of density, interaction strength, and (in the inset) tem- 
perature. It is small deep within the superfluid phase 
and increases to roughly 25% for the largest values of 
U and T shown. Neglecting cubic and quartic terms 
within the Bogoliubov theory relies on the assumption 
of small depletion, and so the conclusions presented here 
are only qualitatively applicable for larger values of U 
and T. (The order of magnitude is consistent with the 
spin- wave analysis of Duric and Lee^) 

For zero temperature, 71r(£ > 0) = and the only 
fluctuation contribution is from the second term within 
the braces [zero-point fluctuations; compare Eq. (|42|) ]. In 
this case, the integrand diverges as u§ | fc | ~ 1 for small |fc|, 
where vq is a coefficient in the expansion of for small 
k (see Appendix IA 3[) . The integral is therefore finite in 
this case. 
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FIG. 7: (Color online) Condensate depletion for a — | as 
a function of interactions U/t (main figure) and temperature 
T/t (inset), where t is the hopping strength. In the main fig- 
ure, T — and the densities are p = 1 (top curve), 2 (middle), 
and 4 (bottom), while in the inset, p = 1 and U/t = 2. The 
depletion is smallest, and hence the approximation best, deep 
in the superfiuid phase, with weak interactions, high density, 
and low temperature. For T > 0, small-momentum cutoffs of 
fen = 0.1 (solid line) and fen = 0.02 (dashed line), in lattice 
units, have been used to remove the logarithmic divergence of 
the depletion integral. 

For T > 0, the Bose-Einstein distribution function be- 
comes T(c|fc|) _1 for small \k\. This results in a logarith- 
mically divergent integral, an instance of the Mermin- 
Wagner-Hohenberg theorem ; 55 ! 56 which states that, in 
two dimensions for nonzero temperature, the continuous 
phase symmetry cannot be broken. In an infinite two- 
dimensional system, there is no true condensate and so 
the 'depletion' is complete. 

In the presence of an external trapping potential, how- 
ever, nonzero temperature condensation is possible even 
in two dimensions^ This can be captured in a crude way 
by applying a small-momentum cutoff ko on the integral 
over fc, with ko ~ R^e ^ where R e g is the effective radius 
of the system in the trap (in units of the lattice spacing). 

If the two-dimensional plane is embedded within a deep 
lattice in the z direction, then hopping in this transverse 
direction can also stabilize the condensate. An appro- 
priate momentum cutoff is then given by ko ~ y/2m*t±, 
where t± is the transverse hopping matrix element, and 
hence the energy scale over which the system appears 
three-dimensional, and m* is the effective mass at the 
minimum of the lowest band in the single-particle dis- 
persion. 

In either case, the depletion integral is finite, with a 
logarithmic dependence on fcp of 

Pio g = ^logfc . (44) 

In the inset of Figure[?l the depletion is shown at nonzero 
temperature, using two different values of ko- The dif- 
ference between the two curves is well approximated for 



most values by Eq. (The exact difference involves 

other terms that are not singular at ko = 0.) To deter- 
mine the cutoff used in the plot, we have assumed that 
the finite system size will be the most important effect, 
and taken 1Z e g ~ 10-50 lattice sites^ 

The depletion calculation also provides a rough esti- 
mate for the boundary of the superfiuid phase, at the 
point where the depletion reaches 100%, although the 
approximation of independent quasiparticles is proba- 
bly not valid at this point. For a = ^, p = 1 and 
T = 0, this gives an estimate of (t/U) c = 0.08, in 
reasonable agreement with the value of (t/U) c = 0.063 
(at the tip of the p = 1 Mott lobe) found using the 
Gutzwiller ansatz . 26 ' 38 ' 39 It should be noted that the lat- 
ter approach, which neglects fluctuations within the Mott 
insulator, generally underestimates {t/U) c ^ 

V. EXPERIMENTAL PREDICTIONS 

The Bogoliubov theory that we have presented for the 
superfiuid phase provides several concrete predictions for 
experiments, most notably for time-of-flight images and 
Bragg spectroscopy. 



A. Time-of-flight images 

As noted above in Section [TTTT the enlarged unit cell of 
the condensate has important consequences for time-of- 
flight images. The corresponding reduction of the Bril- 
louin zone leads to additional Bragg peaks that give a 
clear indication of the formation of spatial order in the 
condensate. The intensity away from these peaks is de- 
termined by bosons excited to states with [fc]<g N ^ by 
thermal and quantum fluctuations. 

In a time-of-flight measurement, the trapping poten- 
tial confining the atoms within the lattice is suddenly 
switched off, causing a rapid expansion. After a fixed 
period of the time, the density profile of the cloud is 
determined, for example by illuminating the atoms and 
measuring the transmitted intensity. If the interactions 
between the atoms during the expansion are sufficiently 
weak, then it can be treated as ballistic, and we will as- 
sume that this is the case throughout. In the absence 
of a magnetic field, the density profile after a fixed time 
of flight measures the original momentum distribution in 
the trap<££ The same is true with a field, apart from some 
modifications that we discuss in the following. 

The time-of-flight images depend on certain details of 
the experiment and, in particular, the means used to 
produce the effective magnetic field. In the case of a 
rotating system, the momentum in the stationary (labo- 
ratory) frame is equal, up to a possible global rotation, 
to the symmetric-gauge canonical momentum in the ro- 
tating frame. 

With a Raman-induced gauge field, the results depend 
on whether the Raman beams remain after release; we 
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assume that they are suddenly switched off simultane- 
ously with the trap, such as in the experiments of Lin et 
al.— The trajectory of an atom is determined by the mo- 
mentum immediately after the gauge field is switched off, 
which is equal, using the sudden approximation, to the 
Landau-gauge canonical momentum before switch-off.»2A 
Within the approximation of a ballistic expansion, the 
time-of-flight image shows the continuum momentum oc- 
cupation N(k), defined by 



N(k) = (¥(k)V(k)) 



(45) 



where ^(fc) is the (continuum) momentum-space annihi- 
lation operator^ The real-space operator ^(r) can, af- 
ter projection to the lowest Bloch band, be expressed in 
terms of the lattice operator bj using the Wannier func- 
tion Wj(r), 



*(r)=]T^(r)&, 



(46) 



In the presence of a magnetic field, this expression cannot 
generally be written as a convolution. 

Combing Eqs. (|35|) and and using Eq. (fTTj) to ex- 
press bj in terms of bk gives 



N(k) 



[ F*{kM){blb k .)F{kM) (47) 



(where, for brevity, the standard integration measure for 
both integrals has been omitted) . The kernel of this dou- 
ble integral transform, analogous to a matrix similarity 
transformation, is given by 



F{k,k') =^ e ife '^ 



d 3 re- ifc - r Wj(r), 



(48) 



where r is integrated over all space. (The Wannier func- 
tion Wj restricts the integral to the neighborhood of the 
two-dimensional plane.) 

The kernel F depends on experimental details, includ- 
ing the optical lattice parameters and the effective gauge 
potential, while theoretical analysis based on the Bosc- 
Hubbard model leads to predictions for the correlation 
function (b k b ko ). We will first outline the form of F 
appropriate to experiments using rotation and Raman- 
induced gauge fields, before giving our results for the 
correlation function based on the Bogoliubov theory. 

In the absence of a magnetic field, the Wannier func- 
tion at site j is a function only of r — Xj , and so can be 

written in the form wj°\r) = u>( )(i — Xj). Shifting the 



integration variable r in Eq. 
evaluated, giving 



allows the sum to be 



F^(fc,fc') = (2^) 2 ,5 2 ([fc-fc']< BL )u; (0) (fc), (49) 

where is the Fourier transform of . This leads 
to the simple result 



N (0 \k) = \w (0 Hk)\ 2 (b 



S b 



(50) 



so the time-of-flight image gives the lattice-momentum 
distribution, with an overall envelope given by the Wan- 
nier function^ 

With nonzero gauge potential ^4, one can instead ex- 
press the Wannier function a o 20 i 22 



Wj{r) = w(i — Xj) exp 



i f dr'- A(r') 



(51) 



where the integral is taken along a straight-line path, as 
in Eq. In this case, the kernel F is no longer diagonal 
in k and k' , and will depend on the appropriate choice 
of gauge. 

In the Landau gauge (with A. parallel to y), one can 
write Ai,(r) = B x r x x, and the kernel is given by 



F L (k,k r ) = J d 3 rc~ ik r w{r)c^ B ^ r y 



x {2ir) 2 5 2 ([k- k' + B x ryy}^) . (52) 
In the symmetric gauge, As(r) — —\B x r, leading to 



Fs(k,k') 



d 3 re-' lkr w(r) 



x (2]r)¥([fe-fc' + ^xr] SL ). (53) 



It should be noted that, in both cases, w(r) differs from 
the Wannier function w^°'(r) in the absence of a mag- 
netic field. 

For our purposes, it is sufficient to note that the ker- 
nels both give contributions to N(k) from a range of 
lattice momenta [fc]<8 L + Sk. The scale is determined by 
\Sk\ < \B\d = 2irad, where d is the characteristic size of 
the Wannier function w. This point-spreading effect can 
be understood as resulting from the position-dependent 
impulse «4.(r) imparted to the atoms when the gauge po- 
tential is switched off. 

We now discuss the form of the correlation function 
(f>fei b ko ), which also depends on the gauge. As in previous 
sections, we will focus on the Landau gauge, appropriate 
for experiments with Raman-induced gauge potentials. 

The symmetry under T£ and implies that {b\ b k ) 
vanishes unless [fei — fc2]s N = 0. We therefore define 

(27r) 2 <5 2 (fe - k')N n ^(k) = (b{ +nX+lY b k , +n , x+tY ) , 

(54) 

for fc, k' € 25n- (The formally infinite factor when k = k' 
corresponds physically to system volume.) 

The dominant contribution to N n £^ n i£i (k) is a delta- 
function peak at k = 0, coming from the first term in 
Eq. HD, 

Nnt.n'Hk) = (2ir) 2 5(k)J2A* ei r in (k)A e , Y ^, n ,(k). 

77' 

(55) 

The corresponding peaks in (b k b k ), at momenta such 
that [fei]<8 N ~ [^2]<8n = 0, are separated by multiples of 
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X and Y. If they are to be resolved in time-of-flight im- 
ages, we require that their separation, 2n/q, be greater 
than the point-spread 2-nad of the kernel F. The con- 
dition is then simply that the Wannier function w(r) be 
well localized compared to the lattice spacing. 

If this condition is satisfied, then the time-of-flight in- 
tensity N(k) consists of sharp peaks near each of the 
reciprocal lattice vectors nX + £Y of the reduced Bril- 
louin zone *Bn- These extra Bragg peaks in fact result 
from two separate physical effects. 

The first is the enlargement of the unit cell of the 
Hamiltonian to q x 1 sites, as a result of the phases 
appearing the hopping term H t . States with momen- 
tum differing by X are therefore mixed at the single- 
particle level, leading to additional Bragg peaks at mo- 
menta 2irxn/q even in the absence of interactions^ The 
observation of such peaks in an experiment is a clear sign 
that the flux per plaquette is at (or sufficiently close to) 
a rational value. 

The second effect occurs only in the presence of inter- 
actions and is due to the spontaneous breaking of spatial 
symmetry in the superfluid. As discussed in Section IlIIl 
the condensate contains contributions from the q degen- 
erate minima of the single-particle dispersion, and there- 
fore enlarges the unit cell to qx q sites. Peaks at momenta 
2?:y£/q are clear indications of the formation of such an 
ordered state. 

Importantly, the kernel for the Landau gauge, 
F^(k,k') in Eq. (|52p . does not change the y component 
of the momentum, and so the point-spreading effect is 
entirely in the x direction. The second class of Bragg 
peaks, resulting from interaction effects, are therefore not 
affected, increasing the likelihood that they can be ob- 
served in experiment. Note that this separation does not 
apply in the symmetric gauge, appropriate for the case of 
rotation, and furthermore that the spacing of the Bragg 
peaks is reduced, as a result of the 2q x 2q unit cell of H t , 
making their observation considerably more challenging. 

Finally, it should be noted that in many cases, includ- 
ing a = 2 and i, the condensate has equal amplitude 
(but not phase) for all values of £. Bragg peaks with 
the same value of n therefore have the same intensity, 
apart from the envelope coming from the on-site Wannier 
wavefunction. This is in contrast to the case of strictly 
vanishing interactions, when any distribution of particles 
between the q minima of the single-particle dispersion 
has equal probability. 

B. Spectroscopy 

Developments in spectroscopic measurements for 
ultracold atomic systems^— have allowed experi- 
mental access to dynamic correlation functions within 
these systems. We consider two such techniques, 
Bragg spectroscop y 62 ' 63 ' 69 and lattice-modulation 



spectroscopy ) 66 ' 6 and describe the information regard- 
ing the quasiparticle spectrum that can be determined 
from both. 

Bragg spectroscop y 62 ' 63 ' 69 involves applying a weak pe- 
riodic perturbation of the form cos(K ■ x — Qt) to the sys- 
tem, using two laser beams at an angle and with frequen- 
cies differing by Q. One then measures, usually through 
time-of-flight imaging, the total momentum or energy im- 
parted to the system. The response is given by the dy- 
namical structure factor S(K, f2), the density correlation 
function in momentum and frequency space. 

Lattice-modulation spectroscop y 66 ' 67 involves oscillat- 
ing the lattice depth at frequency fl; using time-of-flight 
imaging to determine the imparted energy then gives 
S(0, il). It is also possible to measure S(K, fi) at cer- 
tain high-symmetry points K in the lattice Brillouin zone 
2$l by the application of lattices with enlarged periods. 
(Each point in 23l corresponds to a point in Q3n in a 
way that depends on q.) Even if only the point K = 
is accessible, the presence of multiple Hofstadter bands 
should be clear, and the splitting of the Goldstone mode 
from the rest of the first band is also measurable. 

In either case, the coupling to the perturbation can 
be expressed in terms of the momentum-space density 
operator, 

p{k)=[ 7^L^ fe = E e " i!C ^> ( 56 ) 

which, being a function only of rij, is gauge- invariant. 
The dynamic structure factor is given in spectral repre- 
sentation by 

s(K,n) = i J2 c- E * l/T s(E* 2 - E 9l - SI) 

x K^lp^)!^)! 2 , (57) 

where Z = X)* c ~ Eil 'l T is the partition function, and 
|^i 2) are eigenstates of the Hamiltonian T-L with energy 

While S(K, SI) is given by a four-point correlation 
function, it can be factorized into two-point functions 
within the quadratic Bogoliubov theory. In the con- 
densed phase, and assuming depletion is not too large, 
the dominant contribution to the integral in Eq. (|56[) in 
fact comes from the points where either [fc]*8 N = or 
[k — K]<$ N = 0. (For [-K"]<8 N = these cases coincide, 
and there is an extra term in p(K) which, however, con- 
tributes only at ft = 0.) The structure factor is therefore 
given by a two-point correlation function multiplied by 
the condensate density. 

Within this approximation, the density operator can 
be expanded in terms of the operators d kc and d) kC : 
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p([K + NX + LY^J = J2{ri L (K)d K< +u- 2 y<' L lri L (K)]*dl K< } , (58) 



where 

9-1 



r NL 



n,l— 77' 

+ w ^-("+^('+ £ )A [ , +i]a7 ^ [n+JV]4 (0)^ n (-fc)r4,l . (59) 



A subdominant third term involving two operators has been dropped from Eq. 

The dynamic structure factor defined in Eq. (|57p is calculated using the eigenstates of %, which, at the level of 
the quadratic approximation of Section HVl are eigenstates of the occupation numbers of each Bogoliubov mode. The 
matrix elements of Eq. (|58p between any pair of states can be expressed in terms of these occupation numbers, giving 

S([K + NX + LYK.n) = ]T \r c NL (K)\ 2 {s(Q + 6c C )n B (6r C ) + ~ f«rc)[l + MZkc)]} ■ ( 60 ) 



The structure factor at frequency f2 therefore has reso- 
nances at each quasiparticle mode £, allowing the quasi- 
particle spectrum to be measured directly. 

VI. DISCUSSION 

We have studied the effect of a synthetic magnetic field 
on the superfluid phase of bosons in a lattice. Our the- 
oretical approach is based on Bogoliubov theory, which 
determines the condensate configuration and allows inter- 
actions to be taken into account within an expansion in 
terms of fluctuations. We predict broken spatial symme- 
try in the condensed phase, leading to qualitative changes 
compared to the Hofstadter spectrum for noninteracting 
particles. 

This analysis leads to several clear predictions that 
should be testable in experiment. The density modu- 
lations in the superfluid phase, illustrated in Figure 21 
may be directly measurable using recently developed real- 
space imaging techniques! 52 ' 53 Our order-of-magnitude 
estimate for the extent of superfluidity given in Sec- 
tion [IVBl which is in agreement with independent the- 
oretical approaches ) 26 ' 38 ' 39 can be tested in experiments 
analogous to those performed in the absence of a mag- 
netic fields Predictions for spectroscopic measurements 
have been detailed in Section IV Bl 

Our approach also provides predictions for time-of- 
flight imaging, the most well-established technique in 
cold-atom experiments. As described in Section IV A[ we 
predict extra Bragg peaks due to the spatial symmetry 
breaking. In experiments using Raman-induced gauge 
fields, these result from two distinct physical effects. The 
gradient in the applied synthetic vector potential (due 
to a gradient in the physical magnetic field in the ex- 
periments of Lin et ali^i) breaks translation symmetry 



explicitly, leading to an extra set of Bragg peaks in the 
direction of the gradient. By contrast, symmetry un- 
der translation in the perpendicular direction is broken 
spontaneously when the bosons condense, and this leads 
to further Bragg peaks, in the direction of propagation of 
the applied Raman lasers (x in our convention, but y in 
the experiments). The appearance of these latter peaks 
is therefore a clear signature of many-body effects. 

Among the approximations made in the present work is 
the assumption that thermal equilibrium can be reached 
on the time scale of the experiments. Previous studies of 
closely related systems^ have shown that the process of 
vortex formation can exhibit hysteresis, and experiments 
with effective gauge potentials exhibit a considerable de- 
pendence of the vortex density on hold times^ In the 
model considered here, two-body scattering is sufficient 
to populate modes of nonzero £ and hence generate non- 
trivial spatial structures, but further work is required to 
provide quantitative estimates of the rate for these pro- 
cesses. 

We have also neglected the influence of higher lattice 
bands and hopping between pairs of sites other than near- 
est neighbors. Neither is expected to have qualitative 
effects on our conclusions, as long as the magnetic sym- 
metries described in Section III Al are preserved. (This is 
certainly the case with a synthetic magnetic field due to 
rotation or Raman lasers, but not necessarily so when 
hopping phases are induced by other methods^ - — ) As 
already noted in Section II Bl weak interactions between 
bosons on different sites will also have only quantitative 
effects. 

As discussed in Section IIIII the broken spatial sym- 
metry implies the existence of multiple degenerate con- 
figurations in the superfluid phase. This allows for the 
possible formation of real-space domains, especially on 
shorter time scales, upon which the effect of the external 
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trapping potential is likely to be important. 

Besides the simplifications inherent in our starting 
model, our analysis has made the approximation of trun- 
cating the Bogoliubov expansion at quadratic order, 
neglecting interactions between quasiparticles. Conse- 
quences of these interactions include finite quasiparticle 
lifetimes and also the possibility of spectrum termination 
at the point where decay into the two-particle continuum 
is allowed by kinematics. These are likely to have impli- 
cations for spectroscopy experiments; an understanding 
of these is left for future work. 



Acknowledgments 

We thank Ian Spielman, Trey Porto, Chris Foot, Ross 
Williams, and Sarah Al-Assam for helpful discussions. 
This work is supported by JQI-NSF-PFC, ARO-DARPA- 
OLE, and Atomtronics-ARO-MURI. 



Appendix A: Bogoliubov transformation 

In this Appendix, we will show that to diagonalize 
given in Eq. (|37|) . one must find the eigenvalues 
and -vectors of the matrix rjM. k for each k £ Q5n, where 



l,t 7 



Sfp/S. 







77' 



(Al) 



and the matrix product is taken treating both r) and M. k 
as 2q 2 x 2q 2 matrices. While rjMfc is not hermitian, it 
can be shown^ that, since is nonnegative-definite, 
the eigenvalues of rjMfc are all real. Furthermore, for 
k 0, when is positive-definite, the eigenvalues are 
all nonzero and come in pairs with equal magnitude and 
opposite sign. In this case, the q 2 positive eigenvalues 
£kc an d the corresponding eigenvectors V^, 



remains a good symmetry. Instead, define the operator 
V ya = T 2ya 'P representing inversion about the real-space 
point x = yoy (a lattice site if yo is an integer or the 
center of a bond if yo is a half integer). 

Using Eqs. (pH) and (p?6"|) , V yo can be shown to obey 

V yo a kh = e- 2iyok ye- 2i ^cj- 2yoi a { _ k)[ _ e]ql V yQ , (A3) 

so the condensate is invariant under this transformation 
if obeys 



A, 



A 



2i^ UJ -2y £ ^ 



(A4) 



The following assumes that the condensate has an inver- 
sion point, and hence there exists some value of yo for 
which this relation holds. 

Corresponding to this inversion symmetry, and analo- 
gous to the matrix rj, define the 2q 2 x 2q 2 matrix 77, 



71^7,^7' — %+£'],,0^77'e 







It is also convenient^ to define 7, 



lit, 1 



7;*- 7 



Sff'5. 




(A5) 



(A6) 



which has the effect of exchanging creation and annihila- 
tion operators: 70;. = (o^ fc ) T . (The matrices 77, 7r, and 
7 are all both hermitian and unitary, and obey t)tv = tvtj. 

Under the assumption that A( 7 obeys Eq. (|A4[) . 
one can show that 7rMfe7r = M_fc, and furthermore, 
77rM fc 7T7 = M£. These symmetries imply that for ev- 
ery eigenvector of r]M. k with a positive eigenvalue 



£k(, there is a corresponding eigenvector = t^V^, 
with eigenvalue — The corresponding eigenvector of 

r)M_ k is W c _ k = jV C k *, so that £_ hC = £ k( . 



yC _ ( X kt 1 



(A2) 



defined by r;MfeV^ = £k£V£, describe the Bogoliubov 
quasiparticles. 

Section |A 31 treats separately the special case of k = 0, 
allowing us to develop a series expansion for the proper- 
ties near this point. 



1. Inversion symmetry 

The Bogoliubov transformation, which mixes annihila- 
tion operators at momentum k with creation operators 
at — k, requires the existence of an inversion symmetry 
in the condensed phase. Because of broken translational 
symmetry in the presence of a condensate, it is not neces- 
sarily the case that V, inversion about the origin x = 0, 



2. Nonzero momentum 

For nonzero fc, the matrix is positive-definite, and 
all eigenvalues of r/Mfc are real and nonzero. The eigen- 
values then come in pairs of equal magnitude and op- 
posite sign, as claimed previously. It can furthermore be 
shown^i that one can normalize the q 2 vectors so that 



(A7) 



and hence W^' r]Wl = -S cc and w£ rjV k = 0. 

These orthonormality relations immediately lead to 
the results used in Section TlV Al First, they imply that 



K _ 



-t 



rja k 



the operators d k Q 
canonical commutation relations, 



ot_ h rfW''_ h obey the 



d k0 d{, c ,} = (2n) 2 S 2 (k-k')5 cc 



(A8) 
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Second, they lead to the inverse expression 

« fe = E( v ^c + w fc d -*c) ' ( A9 ) 

c 

giving 

«lM fc a fc = J2 ac (d{ C d k C + d k<4<) > (A10) 

C 

from which Eq. P0|) follows. 

3. Near zero momentum 

The invariance of the mean-field energy ho under 
changes of phase of leads to a vanishing eigenvalue 
of the matrix M , which is by assumption the only zero 
eigenvalue (generically the case when no further contin- 
uous symmetries are broken). The corresponding eigen- 
vector is given by 

P^7 = ( ' lAej ) , (All) 

and is obviously also an eigenvector of 77M0 with zero 
eigenvalue. We choose the normalization and phase of P 
so that Ptp = 1 and P = ttP = 7P*. 

It is convenient to define the vector Q satisfying 

r]M Q, = -ivP (A12) 
Q^P = i (A13) 
Q f P = 0, (A14) 

where v is a constant with dimensions of energy. The 
vector Q is specified uniquely by these three equations, 
because M has only one zero eigenvalue, and so its in- 
verse can be defined in the subspace orthogonal to P. 
One can therefore write Q = —wMq 1 t]P, because rfP is 
orthogonal to P, and so is Q by Eq. (|A14|) . Eq. (|A13|) 
simply fixes the normalization of Q, and thus the con- 
stant z/, which can be showr*^ to be positive and given 
by 

v- 1 = P^Mq^P. (A15) 

The vectors P and Q are 'conjugate' in the sense that 
they obey Eq. (|A13|) along with 

Ptr/P = Q f J7Q = 0, (A16) 

so that operators constructed using the vectors P and 
Q obey the commutation relations of momentum and 
position. These operators, however, apply only precisely 
at the zero- measure point fc = 0, and so are not directly 
relevant for the properties in the thermodynamic limit. 

Together with the eigenvectors Vq and Wg corre- 
sponding to nonzero eigenvalues, P and Q span the 2q 2 - 
dimensional vector space, and can therefore be used as 



a basis for a perturbative description of the region near 
fc = 0. For small fc (along x, say), the matrix can 
be expanded as 

Mfc.4 ~ Mo + fc,M w + klM^ , (A17) 

and M^ 1,2 ) can be treated as perturbations. The coef- 
ficient matrices can themselves be calculated using per- 
turbation theory for the wavefunctions i/> 7 „(fc). The re- 
sulting expressions are analytic functions of the momen- 
tum fc, implying the important symmetry property that 
ptMWP = QtM^Q = 0. 

In the following, we restrict to cases where the conden- 
sate is sufficiently symmetric that the phonon velocity is 
isotropic. This requires a symmetry of A under either 1Z 
or I xy (see Appendix [Bj, and as seen in Table U is always 
the case for q < 5. (The extension to general symmetry 
is straightforward.) 

Because of the nonhermitian nature of the matrix 
rjMfc, standard Rayleigh-Schrodinger perturbation the- 
ory cannot be applied directly to this problem, and the 
resulting expressions for the eigenvalues and -vectors are 
not analytic in fc. Instead, the smallest eigenvalue £fci is 
linear in |fc|, 

&i = c|fc| + 0(|fc| 2 ), (A18) 
with phonon velocity c given by 

^^ptM( 2 )P-2V' PtM(1)V °' 2 . (A19) 

As described in Section IA 21 away from fc = all 
eigenvectors have nonzero eigenvalues and can be nor- 
malized so that TjVjj, = -W£ wW^ = 1. As fc 
approaches 0, the vectors V^. and corresponding 
to the smallest eigenvalue both approach P, the unique 
zero-eigenvector of Mo- Since this eigenvector satisfies 
P^P = 0, the normalization of both V^. and must 
diverge as fc — > 0. To leading order, one finds 

Vi =« |fcr 1/2 P + O(|fc| 1/2 ), (A20) 

with coefficient vq = y' vj (2c). The omitted higher- 
order terms maintain the appropriate normalization: 
ptrjVi = i Jo |fc| 1/2 c/V + C(|fc| 3 / 2 ). 

Appendix B: Symmetries and degeneracies 

In Section Hi A[ the magnetic symmetry group (MSG) 
was introduced, and the operators corresponding to var- 
ious elementary operations were defined. In this Ap- 
pendix, we will present in detail the consequences of 
these symmetries for the condensate configurations and 
the quasiparticle spectrum. 

We denote the full group of symmetries of the Hamil- 
tonian W as 0, which includes the translations, rota- 
tions, and reflections considered in Section lll Al As noted 
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in Section IIII1 the mean-field condensate configuration 
breaks a subset of (3; the subgroup of symmetries that 
are preserved will be denoted Sj. Examples are given in 
Table HI which lists the symmetries preserved in the con- 
figurations illustrated in Figure [U We will discuss the 
various properties of the spectrum, including symmetries 
and degeneracies, that result from fj. 

A general spatial transformation, such as an element of 
0, is represented by the (anti)unitary operator 5, under 
which aki-y transforms as 

It should be recalled that reflection operators must be 
combined with time reversal to give symmetries of the 
Hamiltonian, and are hence represented by antiunitary 
operators. The coefficients in Eq. (|B1|) . which can be 
written as a q 2 x q 2 matrix s*,, are diagonal in 7 (with 
possible exceptions at points where two bands touch). 
Preservation of the commutation relations requires that 
Sfc be a unitary matrix, including for antiunitary 5. 

Under the general operation S, the condensate config- 
uration A is mapped to s _1 (0)A(*), with complex con- 
jugation if S is antiunitary. (The inverse matrix arises 
from considering transformations of states versus oper- 
ators.) The mean- field energy ho, defined in Eq. (f33|) . 
is therefore symmetric under the same transformations 
of Agry as the full Hamiltonian is under transformations 
of oo£ 7 - A given configuration that minimizes ho will 
in general have lower symmetry, however, being invari- 
ant only under (a group of mappings isomorphic to) S). 
As argued in Section Mil the superfluid therefore breaks 
spatial symmetries as well as the U(l) phase symmetry. 

As usual, the broken symmetries, comprising the sub- 
set imply the existence of multiple degenerate con- 
figurations. For example, ho is invariant under — > 
Ayt_x\ q ~/ and under A^ — > u~ l Ai 1: corresponding to T x 
and T y respectively [see Eqs. (|2Tj) and (|22|)]. These two 
mappings do not commute, and so there is no config- 
uration A that preserves both symmetries. This leads 
to the conclusion that the number of distinct configu- 
rations that minimize ho is always a multiple of q (by 
the same argument that implies the degeneracy of the 
single-particle states with different £). 

The consequences of the symmetries for the quasi- 
particle spectrum can be determined by considering the 
transformations of the quadratic matrix Mfc. The single- 
particle contribution to given by the first term in 
Eq. (|38p , commutes with any operation S in the full sym- 
metry group © , while the second term results from inter- 
actions with the condensate, and so is symmetric only 
under the elements of Sj. 

These preserved symmetries imply constraints on the 
matrix Mfc. In particular, defining s as in Eq. (|Blj) 
and taking S to be (anti) unitary, one can show that if 



s^AW = crA, then S fc M^ = M fc S fc , where 

s * - (*o l 4) ■ (B2) 

In the presence of a condensate that is symmetric under a 
transformation <S, the quasiparticle energies are therefore 
equal at fc and fc (including in the case where S is antiu- 
nitary). For antiunitary operators, it is convenient to use 
the notation X for the operation of complex conjugation 
followed by multiplication by the matrix S. 

Applied to momenta near k = 0, this leads to con- 
straints on the phonon speed c, calculated in Section TA 31 
In all four cases listed in Table U the symmetry is suf- 
ficient to have isotropic phonon speed, as assumed in 
Eq. (|A18[) . For a = i and ^, the relevant symmetry is 
the (antiunitary) reflection X xy . while for a = j and ^, 
it is the rotation 1Z. (For a = i, there is also symmetry 
under rotation by ^ about a plaquette center, T X TZ.) 

In all cases, the spectrum is of course symmetric only 
under, at most, the fourfold rotation symmetry of the 
square lattice. The symmetry under continuous rotations 
of fc in Eq. (|A18[) is a simple example of an emergent low- 
energy symmetry. 

In many cases (for example a = ^, ^, and ^; see Fig- 
ure @] and Table fl}, t ne condensate configuration pre- 
serves a nontrivial translation symmetry. In particular, 
suppose translation Tr by a displacement R is unbroken, 
where R is not a lattice vector of the enlarged unit cell. 
(Because Tx and T y do not commute, one must specify 
the path to define Tr precisely.) Since translation does 
not change fc, one can define an operator using Eq. (|B2j) 
that commutes with For each fc, the modes ( can 

therefore be labeled according to their eigenvalue under 
Sfc. Modes with different eigenvalues are allowed to have 
(unavoided) crossings, as visible for example in Figure |6] 

1. High-symmetry points 

Points in momentum space separated by the reciprocal 
lattice vectors 2irx and 2iry are physically equivalent, so 
the full Brillouin zone Q5l has the topology of a torus. 
The same applies, with some modifications, to the doubly 
reduced Brillouin zone *Bn- 

With momentum shift operators defined by K^fc = [fc+ 
^-x]<£ L and K y k = [fc + ^-y]<s L , the matrix H nn i(k), 
defined in Eq. (fT7)) . obeys 

H nn i(K x k) = H[n+p] q [n' +p] q {k) (B3) 
H nn ,(K y k) = uj-^ n - n ">H nn ,(k) , (B4) 

so one can extend the definition of -0 7 „(fc) beyond 23n 
by 

Vv^fc) = ^ [n+p]q (k)e ie *^ (B5 ) 
VvOK a fc) = u- pn 4, in {k)c iev ^ k) . (B6) 
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The phases 9y v (k), corresponding to the flux threaded 
through the holes of the torus, are arbitrary apart from 
constraints due to symmetry. These can be found by con- 
sidering the commutation relations of the operators K. x<y 
with each other and with the symmetry operators, and 
by requiring that ip in (fc) be continuous (apart from pos- 
sibly at degeneracy points). The definitions in Eqs. (|B5j) 
and (|B6|) are particularly useful on the boundary of Q3n, 
where the number of constraints on 9^ ,y (k) is larger, and 
especially at the high-symmetry points at the corner (X) 
and edge-center (M) of Q3n- 

The relations between the eigenvectors at k and 
K^k lead to corresponding relations for given by 

K£M fc = M Kf4fe K£, where 



Kf », ,{k) = Sii'Syyuj' 



-«>*(fe) o 



i-yj'Y \ K J - vwo-n"" i Q e i^(fc) 



(B7) 



— ^77' 







e',[e-p] q . 



(B8) 



Appendix C: Analytics for a = \ 

The simplest nontrivial case is a = |, and it is then 
possible to perform many of the calculations analytically. 
In this case, the matrix H(fc) defined in Section Til Bl is 
given by 



. . / , j~\ ,-. , / cos A^7> cos kn 
H(k x x + k y y) = -2t\ x f 
V cos K y — cos k x 



with eigenvalues 



(CI) 



(C2) 



Note that the spectrum has a Dirac cone at the corner of 
25n, where \k x \ = \k y \ = f . At k = 0, the eigenvectors 
are, for bands 7 6 {1,2}, 



6^ = ±2iW cos 2 k x + cos 2 k y 



1 



4 + 2(-l)^^2 



l + (-l)V^ 



(C3) 



These immediately imply that the mode energies are 
equal at points on opposite sides of *Bn- 



2. Kramers degeneracy 

For a = I , there is a twofold degeneracy at every point 
along the line from M to X, as can be seen in Figure [SJ 
and by symmetry at every point on the edge of Q5n- (The 
same degeneracy also occurs for a = j.) 

This is in fact a Kramers degeneracy, and is a conse- 
quence of the symmetry under the glide reflection 7~ x I y 
(for a = 4; the corresponding symmetry for a = j is 
Ty^TxXy). Its action in momentum space is to map k to 
l x k, so that a point k = + k y y, on the line from M 
to X, is mapped to — + k y y, on the opposite side of 
*8n- The quadratic matrix at such a point therefore 
obeys [r^rjMfe] = 0, where 



r fc = Kf^sp" 



(B9) 



with S^ Iv defined in Eq. dB2j) and K£ in Eq. ([B71 
The antiunitary operator has the property 



(BIO) 



which implies, by Kramers' theorem, that all eigenval- 
ues of 'q'M.k are twofold degenerate. Note that, in the 
case a = ■j, this Kramers degeneracy exists despite the 
explicit breaking of time-reversal symmetry by the ap- 
plied magnetic field. (For a = ^, uj = — lis real, and the 
Hamiltonian preserves time-reversal symmetry. The con- 
densate configuration nonetheless breaks this symmetry, 
as shown in Figure 2J) 



The mean-field condensate configuration can be deter- 
mined by calculating ho , given in Eq. (|33[) , and minimiz- 
ing with respect to A^ f at fixed average density. In this 
case, however, the appropriate configuration is more eas- 
ily found by inspection. With the choice Ae t i = i e y/p/2 
and A1.2 = 0, direct calculation shows that the real-space 
wavefunction, given by Eq. (|34p . has uniform magnitude, 



(&,■) = VP exp [(-!)* 



5ir. 
~8~- 



(C4) 



This configuration, in which the condensate is restricted 
to the lower band, therefore has uniform density of p 
particles per lattice site. Calculation of the currents using 
Eq. ([33)1 gives configurations as illustrated in Figure [U 
with the current on each link having magnitude | [J') \ = 
\/2tp. Replacing by its complex conjugate gives the 
equivalent configuration with currents reversed on each 
link. 

Since these configurations have uniform density, they 
simultaneously minimize both terms in Eq. (|33[) globally, 
and are therefore global minima of ho, at fixed density 
p. (They minimize the kinetic energy because they con- 
tain contributions only from the degenerate minima of 
the lowest band, and they minimize the potential energy 
because they have uniform density.) The case a — i is 
unique in this regard, with the condensate configuration 
minimizing both terms simultaneously and so insensitive 
to the value of the interaction strength U. For q > 2, the 
density modulations can be reduced by including higher 
bands in the condensate configuration, and the extent to 
which they contribute is determined by U/t. 

These two configurations are in fact the only minima 
(up to a redundant overall phase rotation), and provide 
an example of the general result that there is always a 
discrete set of degenerate minima, whose number is a 



20 



multiple of q. Either of the translation operators T x and 
T y relates one of the two configurations to the other, up 
to an overall phase (using the transformation of spec- 
ified in Appendix [Bj . 

Both configurations are symmetric under T~ V T X , as 
noted in Table HI with the first obeying 

X> r,r ")^',i =L4i,i, (C5) 

where s 7 " 7 ^ = ( ^ ^ is the matrix defined by Eq. (|B1[) 

for the transformation 7^,7^ (at k = and restricted to 
7 = 1). One can therefore construct the matrix '£ 7 ~ yTl 
according to Eq. (|B2|k it has eigenvalues ±1, allowing the 
quasiparticle modes to be labeled as even or odd under 



Ty%. 



1. Quasiparticle dispersion 

To find the dispersion, one must construct the 8x8 
matrix Mfc given in Eq. ([55)1 . Rather than using the 
single-particle basis labeled by fe, I, and 7, it is somewhat 
easier to find analytic results by starting in the basis of k, 
£, and n as in Eq. (fT6| . While the first term in Eq. (|38|) is 
not diagonal in this basis, the second has a considerably 
simpler expression, as a result of the simple form of the 
on-site interaction in momentum space. 

After transforming to the basis of eigenvectors of 
S 7 " 7 ", the matrix M k splits into two 4x4 blocks, 



M 



±) 



/ U p + 2i(\/2 + cos k x ) 
2t cos k y 
±Up/V2 
\ Up/V2 



2t cos k y ±Up/V2 
Up + 2t(\/2- cos k x ) Up/V2 

Up/V2 Up + 2t(y/2 + cos k x ) 



TUp/V2 



2t cos k,. 



Up 



Up/V2 \ 
TUp/V2 
2t cos k y 

2t(V2-cosk x )J 



(C6) 



corresponding to the eigenvalues ±1. It is then straightforward to calculate the eigenvalues of jyMfc as discussed in 
Section [A II They come in pairs of opposite sign, and so their squares are given by the roots of a quadratic equation. 
The quasiparticle energies are finally given by 



8t 2 + 4V2tUr 



4± 



(m 2 + mV2tUp + 2U 2 p 2 )e 2 k ± m 2 U 2 p 2 cos k x cos h 



y ■ 



(C7) 



where the two choices of ± are independent, giving the 
q 2 = 4 modes of the interacting dispersion. Taking the 
first sign as — and the second as + gives the Goldstone 
mode, which has £ k = \k\ \/ V2pUt + 0(\k\ 3 ) near k = 0. 

Using Eq. (|C7[) , one can confirm the twofold degener- 
acy along the line from M to X established in Section [B~2] 
For these points, cos k x = and the second choice of ± 
is redundant. 



2. Amplitude-phase description of gapless mode 

The small-|fc| dispersion of the Goldstone mode can be 
derived by more elegant means if one restricts to long- 
wavelength fluctuations of the condensate configuration. 
Gradual variations of the real-space wavefunction can be 
parametrized by writing 

&i = <&i>e i *<^l + ^ ) (C8) 

where #j and Qj describe deviations in the phase and 
amplitude respectively, and have canonical commutation 
relations Qj] = i5ij. 



We now rewrite the Hamiltonian H in terms of these 
new degrees of freedom. Assuming Qj <C p and that both 
Qj and $j vary only over distances large compared to the 
lattice scale, one can expand to give 

^o + ^B^-^ + lE^---- (C9) 

(ij) i 

Note that the frustration in T~i t implies that the kinetic 
energy of each link is not separately minimized in the 
mean-field configuration. Each link %— >j therefore con- 
tributes a term linear in §i — $j, but their sum vanishes, 
since the mean-field configuration is a minimum of the 
total kinetic energy. 

Writing this equation in terms of the Fourier compo- 
nents of 'dj and Qj, we obtain 

"^' + /S(> w+ i w ') + -' 

(CIO) 

where the integral is restricted to small |fc| by the as- 
sumption of slowly varying fluctuations. This takes the 
form of a harmonic oscillator for each momentum, so the 

dispersion is = 2 J -j=\k\ 2 x y, in agreement with the 
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result given in the previous section. 
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